setwd("C:/Users/Daniel de Kadt/Dropbox (Personal)/Projects/SA_cohort_effects/replication")
library(foreign)
library(ebal)
library(Matching)
library(xtable)
library(ggplot2)
library(mediation)

data = read.dta("sasas_ddk12.dta")
head(data)
###Varying over force
d = data[which(data$force>-2 & data$force<2 & data$race=="african/black" & (data$vote==1|data$vote==0)),]
d$vote = d$vote*100
obs_diff <- mean(na.omit(d$vote[which(d$treat==1)])) - mean(na.omit(d$vote[which(d$treat==0)]))
sims = 10000

ri_treat_vec=c()
diffs = c()
for(i in 1:sims){
    ri_treat_vec = sample(d$treat,length(d$treat),replace=F)
    diffs[i] <- mean(na.omit(d$vote[which(ri_treat_vec==1)])) - mean(na.omit(d$vote[which(ri_treat_vec==0)]))
}

pval = length(diffs[which(diffs>=obs_diff)])/sims
dd <- with(density(diffs, adjust=2),data.frame(x,y))

pdf(file="sa_1_ri_plot_sasas.pdf", height=6, width=6)
  ggplot(data = dd, mapping = aes(x = x, y = y)) +
    geom_line(color="black") +  layer(data = dd, mapping = aes(x=ifelse(x > obs_diff, x, obs_diff), y=y), geom = "area", geom_params=list(fill="blue",alpha=.3))+
    geom_vline(aes(xintercept=obs_diff), color="blue", linetype="dashed")  +
    scale_y_continuous(limits = c(0,max(dd$y)), name="Density") + 
    scale_x_continuous(name="Difference in Means (Vote)") + 
  ggtitle(paste("1 Year Age Window, Black South Africans ", "(N = ", nrow(d), ")", sep="")) +
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/1.5), text2=paste("Observed =",round(obs_diff,3)), parse=TRUE, size=5)) + 
    geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/4), text2=paste("p-value =",round(pval,3))), size=5)  + 
    theme_bw() +
    theme(axis.line = element_line(colour = "black"),
          panel.border = element_blank(),
          panel.background = element_blank(),
          axis.title.x=element_text(size=12),
          axis.title.y=element_text(size=12,angle=90))
dev.off()

#KENYA
data = read.dta("kenya_plot12.dta")
head(data)
data$vote = data$vote*100
###Varying over force
d = data[which(data$force>-2 & data$force<2 & (data$vote==100|data$vote==0)),]
obs_diff <- mean(na.omit(d$vote[which(d$treat==1)])) - mean(na.omit(d$vote[which(d$treat==0)]))
sims = 10000

ri_treat_vec=c()
diffs = c()
for(i in 1:sims){
  ri_treat_vec = sample(d$treat,length(d$treat),replace=F)
  diffs[i] <- mean(na.omit(d$vote[which(ri_treat_vec==1)])) - mean(na.omit(d$vote[which(ri_treat_vec==0)]))
}

pval = length(diffs[which(diffs>=obs_diff)])/sims
dd <- with(density(diffs, adjust=2),data.frame(x,y))

pdf(file="kenya_1_ri_plot_ab.pdf", height=6, width=6)
ggplot(data = dd, mapping = aes(x = x, y = y)) +
  geom_line(color="black") +  layer(data = dd, mapping = aes(x=ifelse(x > obs_diff, x, obs_diff), y=y), geom = "area", geom_params=list(fill="red",alpha=.3))+
  geom_vline(aes(xintercept=obs_diff), color="red", linetype="dashed")  +
  scale_y_continuous(limits = c(0,max(dd$y)), name="Density") + 
  scale_x_continuous(name="Difference in Means (Vote)") + 
  ggtitle(paste("1 Year Age Window, Kenyans ", "(N = ", nrow(d), ")", sep="")) +
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/1.5), text2=paste("Observed =",round(obs_diff,3)), parse=TRUE, size=5)) + 
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/4), text2=paste("p-value =",round(pval,3))), size=5)  + 
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12,angle=90))
dev.off()

#TANZANIA
data = read.dta("tanzania_plot12.dta")
head(data)
data$vote = data$vote*100
###Varying over force
d = data[which(data$force>-2 & data$force<2 & (data$vote==100|data$vote==0)),]
obs_diff <- mean(na.omit(d$vote[which(d$treat==1)])) - mean(na.omit(d$vote[which(d$treat==0)]))
sims = 10000

ri_treat_vec=c()
diffs = c()
for(i in 1:sims){
  ri_treat_vec = sample(d$treat,length(d$treat),replace=F)
  diffs[i] <- mean(na.omit(d$vote[which(ri_treat_vec==1)])) - mean(na.omit(d$vote[which(ri_treat_vec==0)]))
}

pval = length(diffs[which(diffs <= obs_diff)])/sims
dd <- with(density(diffs, adjust=2),data.frame(x,y))

pdf(file="tanzania_1_ri_plot_ab.pdf", height=6, width=6)
ggplot(data = dd, mapping = aes(x = x, y = y)) +
  geom_line(color="black") +  layer(data = dd, mapping = aes(x=ifelse(x < obs_diff, x, obs_diff), y=y), geom = "area", geom_params=list(fill="red",alpha=.3))+
  geom_vline(aes(xintercept=obs_diff), color="red", linetype="dashed")  +
  scale_y_continuous(limits = c(0,max(dd$y)), name="Density") + 
  scale_x_continuous(name="Difference in Means (Vote)") + 
  ggtitle(paste("1 Year Age Window, Tanzanians ", "(N = ", nrow(d), ")", sep="")) +
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/1.5), text2=paste("Observed =",round(obs_diff,3)), parse=TRUE, size=5)) + 
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/4), text2=paste("p-value =",round(pval,3))), size=5)  + 
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12,angle=90))
dev.off()

#GHANA
data = read.dta("ghana_plot12.dta")
head(data)
data$vote = data$vote*100
###Varying over force
d = data[which(data$force>-2 & data$force<2 & (data$vote==100|data$vote==0)),]
obs_diff <- mean(na.omit(d$vote[which(d$treat==1)])) - mean(na.omit(d$vote[which(d$treat==0)]))
sims = 10000

ri_treat_vec=c()
diffs = c()
for(i in 1:sims){
  ri_treat_vec = sample(d$treat,length(d$treat),replace=F)
  diffs[i] <- mean(na.omit(d$vote[which(ri_treat_vec==1)])) - mean(na.omit(d$vote[which(ri_treat_vec==0)]))
}

pval = length(diffs[which(diffs>=obs_diff)])/sims
dd <- with(density(diffs, adjust=2),data.frame(x,y))

pdf(file="ghana_1_ri_plot_ab.pdf", height=6, width=6)
ggplot(data = dd, mapping = aes(x = x, y = y)) +
  geom_line(color="black") +  layer(data = dd, mapping = aes(x=ifelse(x > obs_diff, x, obs_diff), y=y), geom = "area", geom_params=list(fill="red",alpha=.3))+
  geom_vline(aes(xintercept=obs_diff), color="red", linetype="dashed")  +
  scale_y_continuous(limits = c(0,max(dd$y)), name="Density") + 
  scale_x_continuous(name="Difference in Means (Vote)") + 
  ggtitle(paste("1 Year Age Window, Ghanaians ", "(N = ", nrow(d), ")", sep="")) +
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/1.5), text2=paste("Observed =",round(obs_diff,3)), parse=TRUE, size=5)) + 
  geom_text(mapping=aes(x2,y2,label = text2), data=data.frame(x2=obs_diff, y2=(max(dd$y)/4), text2=paste("p-value =",round(pval,3))), size=5)  + 
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.title.x=element_text(size=12),
        axis.title.y=element_text(size=12,angle=90))
dev.off()

